Polyfem Tutorial¶
This is a jupyter notebook. The “real” notebook can be found here.
Polyfem relies on 3 main objects:
Settingsthat contains the main settings such discretization order (e.g., P_1 or P_2), material parameters, formulation, etc.Problemthat describe the problem you want to solve, that is the boundary conditions and right-hand side. There are some predefined problems, such asDrivenCavity, or generic problems, such asGenericTensor.Solverthat is the actual FEM solver.
The usage of specific problems is indented for benchmarking, in general you want to use the GenericTensor for tensor-based PDEs (e.g., elasticity) or GenericScalar for scalar PDEs (e.g., Poisson).
A typical use of Polyfem is:
settings = polyfempy.Settings() # set necessary settings # e.g. settings.discr_order = 2 problem = polyfempy.GenericTensor() # or any other problem # set problem related data # e.g. problem.set_displacement(1, [0, 0], [True, False]) settings.set_problem(problem) #now we can create a solver and solve solver = polyfempy.Solver() solver.settings(settings) solver.load_mesh_from_path(mesh_path) solver.solve()
Note 1: for legacy reasons Polyfem always normalizes the mesh (i.e., rescale it to lay in the [0,1]^d box, you can use setting.normalize_mesh = False to disable this feature.
Note 2: the solution u(x) of a FEM solver are the coefficients u_i you need to multiply the bases \varphi_i(x) with: $$ u(x)=\sum u_i \varphi_i(x). $$ The coefficients u_i are unrelated with the mesh vertices because of reordering of the nodes or high-order bases. For instance P_2 bases have additional nodes on the edges which do not exist in the mesh.
For this reason Polyfem uses a visualization mesh where the solution is sampled at the vertices. This mesh has two advantages: 1. it solves the problem of nodes reordering and additional nodes in the same way 2. it provides a “true” visualization for high order solution by densely sampling each element (a P_2 solution is a piecewise quadratic function which is visualized in a picewise linear fashion, thus the need of a dense element sampling).
To control the resolution of the visualization mesh use settings.vismesh_rel_area.
Examples¶
Some imports for plotting
import plotly.offline as plotly import plotly.graph_objs as go import plotly.figure_factory as ff website_mode = False #Necessary for the notebook plotly.init_notebook_mode(connected=not website_mode)
algebra
import numpy as np
stuff for the animation
import ipywidgets as widgets from ipywidgets import interact
and finallypolyfempy
import polyfempy as pf
Plotting utility¶
This code is not particularly interesting.
It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.
def create_plot_mesh_and_layout(vertices, connectivity, function): x = vertices[:,0] y = vertices[:,1] if vertices.shape[1] == 3: z = vertices[:,2] else: z = np.zeros(x.shape, dtype=x.dtype) if connectivity.shape[1] == 3: f = connectivity else: #Convert a tet-mesh into face based triangles. f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64) for i in range(len(connectivity)): f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]]) f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]]) f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]]) f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]]) mesh = go.Mesh3d(x=x, y=y, z=z, i=f[:,0], j=f[:,1], k=f[:,2], intensity=function, flatshading=function is not None, colorscale='Viridis', contour=dict(show=True, color='#fff'), hoverinfo="all") layout = go.Layout(scene=go.layout.Scene( aspectmode='data', xaxis=dict( autorange=True, showgrid=False, zeroline=False, showline=False, ticks='', showticklabels=False ), yaxis=dict( autorange=True, showgrid=False, zeroline=False, showline=False, ticks='', showticklabels=False ), zaxis=dict( autorange=True, showgrid=False, zeroline=False, showline=False, ticks='', showticklabels=False ) )) return mesh, layout
def plot(vertices, connectivity, function, camera=None): mesh, layout = create_plot_mesh_and_layout(vertices, connectivity, function) if camera is not None: layout.scene.camera = camera fig = go.Figure(data=[mesh], layout=layout) if website_mode: return fig plotly.iplot(fig)
Creates a quad mesh of n_pts x n_pts in the form of a regular grid
def create_quad_mesh(n_pts): extend = np.linspace(0,1,n_pts) x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy') pts = np.column_stack((x.ravel(), y.ravel())) faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32) index = 0 for i in range(n_pts-1): for j in range(n_pts-1): faces[index, :] = np.array([ j + i * n_pts, j+1 + i * n_pts, j+1 + (i+1) * n_pts, j + (i+1) * n_pts ]) index = index + 1 return pts, faces
Plate hole¶
This is the python version of the plate with hole example explained here.
Set the mesh path
mesh_path = "plane_hole.obj"
create settings
settings = pf.Settings()
pick linear P_1 elements. If the mesh would be a quad it would be Q_1
settings.discr_order = 1
normalize the mesh to be in [0,1]^2
settings.normalize_mesh = True
and choose Young’s modulus and poisson ratio
settings.set_material_params("E", 210000) settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
Next we setup the problem
problem = pf.GenericTensor()
sideset 1 has zero displacement in x
problem.set_displacement(1, [0, 0], [True, False])
sideset 4 has zero displacement in y
problem.set_displacement(4, [0, 0], [False, True])
sideset 3 has a force of [100, 0] applied
problem.set_force(3, [100, 0])
fianally set the problem
settings.set_problem(problem)
Solve!
solver = pf.Solver() solver.settings(settings) solver.load_mesh_from_path(mesh_path) solver.solve()
[2019-05-28 18:38:02.409] [polyfem] [info] Loading mesh... [2019-05-28 18:38:02.409] [geogram] [info] Loading file plane_hole.obj... [2019-05-28 18:38:02.419] [geogram] [info] (FP64) nb_v:8549 nb_e:0 nb_f:16797 nb_b:299 tri:1 dim:3 [2019-05-28 18:38:02.419] [geogram] [info] Attributes on vertices: point[3] [2019-05-28 18:38:02.428] [polyfem] [info] mesh bb min [0.500183, 0.5], max [100.5, 50.5] [2019-05-28 18:38:02.428] [polyfem] [info] took 0.0190063s [2019-05-28 18:38:02.430] [polyfem] [info] Building isoparametric basis... [2019-05-28 18:38:02.450] [polyfem] [info] Computing polygonal basis... [2019-05-28 18:38:02.450] [polyfem] [info] took 1.6903e-05s [2019-05-28 18:38:02.451] [polyfem] [info] hmin: 0.420699 [2019-05-28 18:38:02.451] [polyfem] [info] hmax: 1.875 [2019-05-28 18:38:02.451] [polyfem] [info] havg: 0.831949 [2019-05-28 18:38:02.451] [polyfem] [info] took 0.0202445s [2019-05-28 18:38:02.451] [polyfem] [info] flipped elements 0 [2019-05-28 18:38:02.451] [polyfem] [info] h: 1.875 [2019-05-28 18:38:02.451] [polyfem] [info] n bases: 8549 [2019-05-28 18:38:02.451] [polyfem] [info] n pressure bases: 0 [2019-05-28 18:38:02.451] [polyfem] [info] Assigning rhs... [2019-05-28 18:38:02.452] [polyfem] [info] took 0.000734283s [2019-05-28 18:38:02.452] [polyfem] [info] Assembling stiffness mat... [2019-05-28 18:38:02.484] [polyfem] [info] took 0.0320988s [2019-05-28 18:38:02.484] [polyfem] [info] sparsity: 236956/292341604 [2019-05-28 18:38:02.484] [polyfem] [info] Solving LinearElasticity with [2019-05-28 18:38:02.484] [polyfem] [info] Hypre...
Get the solution
[pts, tets, disp] = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
top_camera = dict( up=dict(x=0, y=1, z=0), center=dict(x=0, y=0, z=0), eye=dict(x=0, y=0, z=3.5) ) plot(vertices, tets, mises, top_camera)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use P_1 elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the “real” mises just call
mises = solver.get_sampled_mises() plot(vertices, tets, mises, top_camera)